Three Rivers Parks species lists

Load packages

library(tidyverse)
library(ggplot2)
library(sf)
library(rgbif)
library(here)
library(CoordinateCleaner)
library(sp)
library(DT)

"%!in%" <- function(x, y) !("%in%"(x, y))

Load park boundaries

Load park boundaries from which to compile species occurrence data. This is a manually edited version of the boundary file provided by TRPD which resolves some topology issues. It is also slightly buffered (30 m) to capture records immediately outside the legal park boundaries (e.g. records from adjacent roads).

parks_sf <- readRDS(here("Data", "Wildlife", "eBird", "parks_outlines_ManualEdits_30mBufferDissolve.RDS")) %>%
  st_transform(crs = 4326) %>%
  select(LabelName) %>%
  rename(Park_Name = LabelName)
park_names <- parks_sf$Park_Name

# Not run now, but this code generates the WKT from which a spatial search is performed on the GBIF website, using the convert shape to bounding box function.
# parks_sp <- as_Spatial(parks_sf)
# parks_wkt <- sp_convert(parks_sp) # convert sp object to WKT for rgbif::occ_search geometry argument

Load species occurrence data

The Global Biodiversity Information Facility compiles species occurrence data from many sources, including museums and participatory science platforms like eBird and iNaturalist.

GBIF data queries are saved on the GBIF servers and get a unique identifier and DOI. The code below inputs a GBIF data query for all animals in the rectangular area bounding all Three Rivers parks. It then filters down the data to only records that come from within the park boundaries and pass a few basic quality filters. Information on the specific query can be found on this page. This version uses a query made on 2023-07-19.

# Define where to save the filtered data. If this file already exists the script will simply load it instead of re-downloading and filtering.
filtered_occ_dat_file <- here("Data", "Wildlife", "GBIF", "0101388-230530130749713_filtered.RDS")

if (!file.exists(filtered_occ_dat_file)) {
  occurence_data_all <- occ_download_get("0101388-230530130749713", path = here("Data", "Wildlife", "GBIF"), overwrite = FALSE) %>% occ_download_import()

  occurence_data_parks <- occurence_data_all %>%
    setNames(tolower(names(.)))

  # convert occurrence data to spatial features and filter to park boundaries
  occurence_data_parks <- occurence_data_parks %>%
    st_as_sf(coords = c("decimallongitude", "decimallatitude"), crs = 4326, remove = FALSE)
  occurence_data_parks <- occurence_data_parks %>%
    st_filter(parks_sf)

  # add park name each record falls within to dataframe. Doing this in a loop for each park seems to be faster than st_intersect().
  parks_mp <- parks_sf %>%
    group_by(Park_Name) %>%
    summarise(geometry = st_combine(geometry))

  dat_parks_int <- data.frame()
  for (i in 1:nrow(parks_mp)) {
    park_name <- parks_mp$Park_Name[[i]]
    park <- parks_mp[i, ]
    park_dat <- st_filter(occurence_data_parks, park)
    park_dat$Park_Name <- park_name
    dat_parks_int <- rbind(dat_parks_int, park_dat)
  }

  occurence_data_parks <- dat_parks_int
  rm(dat_parks_int)

  ## Data filtering

  # in previous versions of this analysis I used rgbif's occ_search() to get data, which could be filter ed based on predefined issues. This doesn't seem to work with he version imported by occ_download_import(), but saving code for future reference.
  # Filter pre-defined GBIF issues
  # issues<-gbif_issues()
  # occurence_no_issues <- temp %>%
  #   occ_issues(-ccm, -cdiv, -cdout, -cdrepf, -cdreps, -cucdmis, -preneglat, -preneglon, -preswcd, -rdativ, -rdatm, -rdatunl, -txmathi, -txmatnon, -typstativ, -zerocd, -interr) # May want to add txmatfuz
  # occurence_data_filtered <- occurence_no_issues$data

  # Other filters
  occurence_data_filtered <- occurence_data_parks %>%
    filter(!is.na(year)) %>% # remove records with no year
    filter(!is.na(species)) %>% # remove records with no species
    filter(species != "") %>% # remove records with no species
    filter(occurrencestatus == "PRESENT") %>% # only records indicating presence
    filter(!is.na(decimallongitude)) %>% # drop records without longitude
    filter(!is.na(decimallatitude)) %>% # drop records without latitude
    filter(!basisofrecord %in% c("FOSSIL_SPECIMEN", "LIVING_SPECIMEN")) %>% # drop fossils and "living specimens", which are usually captive things
    filter(year >= 1500) %>% # remove really old records
    # filter(coordinateprecision >= 0.1 | is.na(coordinateprecision)) %>% # remove records with large coordinate uncertainty. (TK NOT WORKING)
    filter(coordinateuncertaintyinmeters < 60000 | is.na(coordinateuncertaintyinmeters)) %>% # remove records with large coordinate uncertainty. Elm Creek is largest park by area listed on Wikipedia at 21.2 km2, which as a circle would have a radius of ~5.2 km. Uncertainty must be less than this. TK make this dependent on size of actual park?
    filter(!coordinateuncertaintyinmeters %in% c(301, 3036, 999, 9999)) %>% # remove records with known default coordinate values
    filter(!decimallatitude == 0 | !decimallongitude == 0) %>% # remove point plotted along the prime meridian or equator
    glimpse() # look at results of pipeline

  # remove records with certain words and phrases in certain fields
  problem_words_locality <- "\\bzoo\\b|aquarium" # remove zoo (whole word) and aquarium from verbatim locality

  occurence_data_filtered <- occurence_data_filtered %>%
    select(-c(datasetkey, stateprovince, publishingorgkey, elevation, elevationaccuracy, depth, depthaccuracy, recordnumber, dateidentified, license, rightsholder, typestatus, establishmentmeans, mediatype))

  saveRDS(occurence_data_filtered, filtered_occ_dat_file)
} else {
  occurence_data_filtered <- readRDS(filtered_occ_dat_file)
}

# apply manual changes to individual records

# eastern blacknose dace should be western backnose dace (split)
occurence_data_filtered$species[which(occurence_data_filtered$species == "Rhinichthys atratulus")] <- "Rhinichthys obtusus"

All parks

At least 988 animal species have been documented in Three Rivers Parks. These species include 295 birds, 38 fishes, 29 mammals, 11 reptiles, and 10 amphibians, and 570 arthropods. The arthropod species include 156 butterflies/moths; 80 bees/wasps/ants/sawflies; 75 dragonflies; 71 beetles; 65 flies; 46 spiders; 36 true bugs; and 21 grasshoppers.

The full list of documented species is provided below. This table can be filtered or searched to find particular taxa.

Individual parks

This version of the table can be filtered to get the species list for any individual park. Or, an individual species can be searched to determine which parks it has been documented in.

All parks

Species possibly extirpated from all parks:

Individual parks

Species possibly extirpated from individual parks (Park Reserves only):

Species on TRPD lists but not in GBIF

The table below lists species of mammals, herpetofauna, and fish appearing on Three Rivers species lists but that are not documented in the GBIF database. Many or most of these “missing” species likely do occur in the parks but have gone undocumented in the databases reporting to GBIF. TRPD Wildlife staff could consider adding records of these species to iNaturalist or directly to GBIF to help fill these data gaps. The public could also be solicited to submit documentation for these species.

# source lists
source(here("Scripts", "976_TRPD_species_lists_as_dataframes.R"))

# function to match species names on TRPD lists with GBIF taxonomy
accepted_name <- function(name) {
  temp <- name_backbone(name)
  temp$species
}

# add GBIF accepted name to TRPD list for mammals, herps, fish
mammals_data$accepted_name <- sapply(mammals_data$scientific_name, accepted_name)
herpetofauna_data$accepted_name <- sapply(herpetofauna_data$scientific_name, accepted_name)
fish_data$accepted_name <- sapply(fish_data$scientific_name, accepted_name)
fish_data$accepted_name <- sapply(fish_data$scientific_name, accepted_name)
hymenoptera_data$accepted_name <- sapply(hymenoptera_data$scientific_name, accepted_name)

# identify species not documented in GBIF within parks (but on TRPD species list)
mammals_not_in_GBIF <- setdiff(mammals_data$accepted_name, species_summary$species)
mammals_not_in_GBIF <- mammals_data[which(mammals_data$accepted_name %in% mammals_not_in_GBIF), ]
mammals_not_in_GBIF %>%
  select(common_name, accepted_name) %>%
  rename(
    "common name" = common_name,
    "scientific name" = accepted_name
  )

herpetofauna_not_in_GBIF <- setdiff(herpetofauna_data$accepted_name, species_summary$species)
herpetofauna_not_in_GBIF <- herpetofauna_data[which(herpetofauna_data$accepted_name %in% herpetofauna_not_in_GBIF), ]
herpetofauna_not_in_GBIF %>%
  select(common_name, accepted_name) %>%
  rename(
    "common name" = common_name,
    "scientific name" = accepted_name
  )

fish_not_in_GBIF <- setdiff(fish_data$accepted_name, species_summary$species)
fish_not_in_GBIF <- fish_data[which(fish_data$accepted_name %in% fish_not_in_GBIF), ]
fish_not_in_GBIF %>%
  select(common_name, accepted_name) %>%
  rename(
    "common name" = common_name,
    "scientific name" = accepted_name
  )

hymenoptera_not_in_GBIF <- setdiff(hymenoptera_data$accepted_name, species_summary$species)
hymenoptera_not_in_GBIF <- hymenoptera_data[which(hymenoptera_data$accepted_name %in% hymenoptera_not_in_GBIF), ]

species_not_on_GBIF <- rbind(mammals_not_in_GBIF, herpetofauna_not_in_GBIF, fish_not_in_GBIF)

In total there are 87 species of mammals, herpetofauna, and fishes that are listed as occurring in the parks on TRPD lists that do not have any records in the GBIF database. Most of the “missing” mammals are small rodents, fossorial species, or nocturnal species, which are difficult to observe/document without dedicated surveys/trapping. Some small-mammal trapping has been conducted in the parks. Data from these efforts could be reviewed and submitted to a GBIF data source to improve the records. Fish are systematically underrepresented in citizen science databases like iNaturalist, which likely explains the large number of “missing” species in that group. The missing herpetofauna includes at least one species with documented iNaturalist records in the park (plains hog-nosed snake, Heterodon platirhinos), but since they are a sensitive species, the records are obscured and are not included in the public data download.

Species in GBIF but not on TRPD list

In this version of the table we display species in the GBIF database that are not included on TRPD species lists maintained by the Wildlife Management Department.

There are only 2 species of mammals, fish, or herpetofauna that have been documented in the parks that do not already appear on the TRPD species lists:

• Domestic cat (Felis catus): an introduced species. Since the mammal list maintained by TRPD contains other introduced species, it could be appropriate to also add cat. The single record of cat in this database comes from the southern edge of Murphy-Hanrehan Park.

• Ouachita Map Turtle (Graptemys ouachitensis): represented by a single record from Medicine Lake in French Regional Park. Almost certainly a released pet and not known to be part of a self-sustaining population.

This exercise is perhaps more useful for taxonomic groups where park species lists are more uncertain or derived from general county-level datasets. For example there are approximately 38 species in the order Hymenoptera (bees/wasps/sawflies) that are not currently on the TRPD list but that have been documented in the parks on GBIF.

Species accumulation curves

Recording of new species increased rapidly in late 2010s–perhaps starting to level off in some parks.

# Count new species in each year by park
species_accumulation_indvidual_parks <- species_list_parks %>%
  filter(str_detect(Park_Name, "Reserve")) %>%
  group_by(Park_Name, min_year) %>%
  summarise(n = n()) %>%
  arrange(min_year) %>%
  mutate(running_sp_total = cumsum(n))

# Determine the maximum value for each park for setting legend
max_values <- species_accumulation_indvidual_parks %>%
  group_by(Park_Name) %>%
  summarise(max_value = max(running_sp_total))

# Sort the Park_Name factor levels by their maximum values
species_accumulation_indvidual_parks$Park_Name <- factor(species_accumulation_indvidual_parks$Park_Name, levels = max_values$Park_Name[order(max_values$max_value)])

# Create plot
ggplot(species_accumulation_indvidual_parks, aes(x = min_year, y = running_sp_total, group = Park_Name, color = Park_Name)) +
  geom_point() +
  geom_line() +
  ggthemes::theme_clean() +
  labs(
    y = "Cumulative species count",
    x = "Year"
  ) +
  theme(legend.position = "right") +
  guides(color = guide_legend(reverse = TRUE)) +
  coord_cartesian(xlim = c(1980, 2023))

Documented species richness by park

Plots of total number of species by park

sp_richness_plot <- species_list_parks %>%
  mutate(group = case_when(
    phylum %in% c("Annelida", "Mollusca", "Rotifera") ~ "Other invertebrates",
    phylum == "Arthropoda" & class %!in% c("Insecta", "Arachnida") ~ "Other invertebrates",
    class == "Arachnida" ~ "Arachnids",
    class == "Insecta" & order %!in% c("Lepidoptera", "Odonata", "Hymenoptera", "Diptera", "Coleoptera") ~ "Other invertebrates",
    class == "Insecta" & order %in% c("Lepidoptera", "Odonata", "Hymenoptera", "Diptera", "Coleoptera") ~ order,
    class %in% c("Squamata", "Testudines", "Amphibia") ~ "Herpetofauna",
    TRUE ~ as.character(class)
  )) %>%
  mutate(group = as.factor(group)) %>%
  mutate(group = fct_relevel(group, "Arachnids", "Diptera", "Coleoptera", "Odonata", "Hymenoptera", "Lepidoptera", "Other invertebrates", "Mammalia", "Actinopterygii", "Herpetofauna", "Aves")) %>%
  # filter(str_detect(Park_Name, "Reserve")) %>%
  group_by(Park_Name, group) %>%
  summarise(n = n()) %>%
  ggplot() +
  geom_bar(aes(x = reorder(Park_Name, n, sum), y = n, fill = group), position = "stack", stat = "identity") +
  ggthemes::theme_clean() +
  labs(
    y = "Documented species richness",
    x = ""
  ) +
  theme(
    legend.position = "right",
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)
  ) +
  scale_x_discrete(limits = rev)

sp_richness_plot

Version of same plot as above, but with an indication of sampling effort. Black dots show an index of total number of bird records. Grey dots are total number of non-bird records.

max_sp_park <- species_list_parks %>%
  st_drop_geometry() %>%
  group_by(Park_Name) %>%
  summarize(n = n()) %>%
  pull(n) %>%
  max()

record_count <- occurence_data_filtered %>%
  st_drop_geometry() %>%
  group_by(Park_Name) %>%
  summarise(n = n())

record_count_nobirds <- occurence_data_filtered %>%
  filter(class != "Aves") %>%
  st_drop_geometry() %>%
  group_by(Park_Name) %>%
  summarise(n = n())

record_count_birds <- occurence_data_filtered %>%
  filter(class == "Aves") %>%
  st_drop_geometry() %>%
  group_by(Park_Name) %>%
  summarise(n = n())

sp_richness_plot +
  geom_point(data = record_count_birds, aes(x = Park_Name, y = n / (max(record_count_birds$n) / max_sp_park))) +
  geom_point(data = record_count_nobirds, aes(x = Park_Name, y = n / (max(record_count_nobirds$n) / max_sp_park)), color = "grey")

This shows that Crow-Hassan has high species count despite relatively low sampling effort. Murphy way more eBird use than iNat use…

sessionInfo()

## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8  LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] CoordinateCleaner_3.0.1 rgbif_3.7.8             MASS_7.3-58.1          
##  [4] viridis_0.6.4           viridisLite_0.4.2       effects_4.2-2          
##  [7] carData_3.0-5           knitr_1.45              snakecase_0.11.1       
## [10] DHARMa_0.4.6            glmmTMB_1.1.8           performance_0.10.8     
## [13] insight_0.19.7          gridExtra_2.3           jagsUI_1.5.2           
## [16] png_0.1-8               transformr_0.1.3        gifski_1.12.0-2        
## [19] gganimate_1.0.8         ggspatial_1.1.9         ggrepel_0.9.4          
## [22] RVAideMemoire_0.9-83-7  pairwiseAdonis_0.4.1    cluster_2.1.4          
## [25] goeveg_0.6.5            vegan_2.6-4             lattice_0.20-45        
## [28] permute_0.9-7           mapview_2.11.2          ggtext_0.1.2           
## [31] ratelimitr_0.4.1        rvest_1.0.3             trelliscopejs_0.2.6    
## [34] plotly_4.10.3           auk_0.7.0               readxl_1.4.3           
## [37] kableExtra_1.3.4        ggthemes_5.0.0          forcats_1.0.0          
## [40] stringr_1.5.1           purrr_1.0.2             readr_2.1.4            
## [43] tidyr_1.3.0             tibble_3.2.1            ggplot2_3.4.4          
## [46] tidyverse_2.0.0         lubridate_1.9.3         dplyr_1.1.4            
## [49] DT_0.31                 sf_1.0-14               rgdal_1.6-4            
## [52] sp_2.1-2                RODBC_1.3-23            here_1.0.1             
## 
## loaded via a namespace (and not attached):
##   [1] estimability_1.4.1      coda_0.19-4             ragg_1.2.6             
##   [4] bit64_4.0.5             multcomp_1.4-25         data.table_1.14.8      
##   [7] rpart_4.1.19            doParallel_1.0.17       generics_0.1.3         
##  [10] leaflet_2.2.1           terra_1.7-55            cowplot_1.1.1          
##  [13] TH.data_1.1-2           commonmark_1.9.0        proxy_0.4-27           
##  [16] bit_4.0.5               tzdb_0.4.0              webshot_0.5.5          
##  [19] xml2_1.3.5              httpuv_1.6.12           wk_0.9.1               
##  [22] assertthat_0.2.1        oai_0.4.0               xfun_0.41              
##  [25] hms_1.1.3               jquerylib_0.1.4         satellite_1.0.4        
##  [28] evaluate_0.23           promises_1.2.1          fansi_1.0.5            
##  [31] progress_1.2.3          DBI_1.1.3               htmlwidgets_1.6.4      
##  [34] stats4_4.2.2            ellipsis_0.3.2          crosstalk_1.2.1        
##  [37] backports_1.4.1         survey_4.2-1            markdown_1.12          
##  [40] epuRate_0.1             vctrs_0.6.5             geosphere_1.5-18       
##  [43] rnaturalearth_0.3.4     cachem_1.0.8            withr_2.5.2            
##  [46] triebeard_0.4.1         ggh4x_0.2.6             checkmate_2.3.0        
##  [49] vroom_1.6.4             emmeans_1.8.9           prettyunits_1.2.0      
##  [52] mclust_6.0.1            svglite_2.1.2           dotCall64_1.1-1        
##  [55] lazyeval_0.2.2          crayon_1.5.2            leaflet.providers_2.0.0
##  [58] crul_1.4.0              pkgconfig_2.0.3         labeling_0.4.3         
##  [61] units_0.8-5             tweenr_2.0.2            nlme_3.1-160           
##  [64] nnet_7.3-18             rlang_1.1.2             lifecycle_1.0.4        
##  [67] sandwich_3.0-2          httpcode_0.3.0          cellranger_1.1.0       
##  [70] rprojroot_2.0.4         urltools_1.7.3          Matrix_1.6-4           
##  [73] raster_3.6-26           boot_1.3-28             zoo_1.8-12             
##  [76] base64enc_0.1-3         whisker_0.4.1           gap.datasets_0.0.6     
##  [79] KernSmooth_2.23-20      spam_2.10-0             classInt_0.4-10        
##  [82] s2_1.1.4                brew_1.0-8              scales_1.3.0           
##  [85] lpSolve_5.6.19          magrittr_2.0.3          plyr_1.8.9             
##  [88] compiler_4.2.2          lme4_1.1-35.1           cli_3.6.1              
##  [91] pbapply_1.7-2           TMB_1.9.9               htmlTable_2.4.2        
##  [94] Formula_1.2-5           mgcv_1.8-41             tidyselect_1.2.0       
##  [97] stringi_1.8.2           textshaping_0.3.7       DistributionUtils_0.6-1
## [100] highr_0.10              mitools_2.4             yaml_2.3.7             
## [103] rebird_1.3.0            grid_4.2.2              sass_0.4.7             
## [106] tools_4.2.2             timechange_0.2.0        rstudioapi_0.15.0      
## [109] uuid_1.1-1              foreach_1.5.2           foreign_0.8-83         
## [112] rjags_4-15              leafpop_0.1.0           farver_2.1.1           
## [115] digest_0.6.33           shiny_1.8.0             qgam_1.3.4             
## [118] autocogs_0.1.4          Rcpp_1.0.11             MCMCvis_0.16.3         
## [121] gridtext_0.1.5          later_1.3.1             httr_1.4.7             
## [124] Rdpack_2.6              colorspace_2.1-0        splines_4.2.2          
## [127] fields_15.2             systemfonts_1.0.5       xtable_1.8-4           
## [130] jsonlite_1.8.7          nloptr_2.0.3            leafem_0.2.3           
## [133] gap_1.5-3               R6_2.5.1                Hmisc_5.1-1            
## [136] pillar_1.9.0            htmltools_0.5.7         mime_0.12              
## [139] glue_1.6.2              fastmap_1.1.1           minqa_1.2.6            
## [142] class_7.3-20            codetools_0.2-18        maps_3.4.1.1           
## [145] mvtnorm_1.2-4           utf8_1.2.4              bslib_0.6.1            
## [148] numDeriv_2016.8-1.1     curl_5.1.0              unmarked_1.3.2         
## [151] survival_3.4-0          rmarkdown_2.25          munsell_0.5.0          
## [154] e1071_1.7-13            iterators_1.0.14        gtable_0.3.4           
## [157] rbibutils_2.2.16